\documentclass[a4paper]{article}
\usepackage{amsmath,amssymb,amsfonts,latexsym,graphicx}
\usepackage[boldfont,slantfont]{xeCJK}
\usepackage{xcolor}
\usepackage{enumerate}
\setCJKmainfont{AR PL UMing CN}


\parindent=0mm
\textheight=20cm \textwidth=14cm
\parskip=0mm
\renewcommand\baselinestretch{1.2}
\oddsidemargin30pt \evensidemargin30pt


\begin{document}

\baselineskip 18pt
\parskip 10pt
\parindent 0em

\def\chntoday{\the\year~年~\the\month~月~\the\day~日}

\renewcommand\figurename{图}\renewcommand\tablename{表}
\renewcommand{\refname}{参考文献}
\renewcommand\abstractname{\bf\large 摘~~要}

\title{计算流体力学课程上机作业(二)}
\author{叶时炜}
\date{\chntoday}
\maketitle

\begin{abstract}
一维线性Euler方程组的初边值问题。用Lax-Friedrichs格式和MacCormack格式求不同初值条件下的程序。
\end{abstract}

\section{问题}

已知一维Euler方程组：
\begin{displaymath}\label{eq:euler}
  \left\{
  \begin{array}{c}
    \left(
    \begin{array}{c}
      \rho\\
      \rho u\\
      E\\
    \end{array}
    \right)_t + 
    \left(
    \begin{array}{c}
      \rho u\\
      \rho u^2 +p\\
      u(E+p)\\
    \end{array}
    \right)_x=0,\\
    p=(\gamma-1)(E-\frac{1}{2}\rho u^2),\gamma=1.4\\
  \end{array}
  \right.
\end{displaymath}
用Lax-Friedrichs格式和MacCormack格式求下列初边值问题。
\begin{description}
\item[问题3]初始条件：
  \begin{displaymath}
    (\rho, u, p)(x,0)=\left\{
    \begin{array}{ll}
      (1, 0, 10^3),&0\leq x <0.1,\\
      (1, 0, 10^{-2}),&0.1\leq x <0.9,\\
      (1, 0, 10^2),&0.9\leq x <1,\\
    \end{array}
    \right.
  \end{displaymath}
  计算区间为[0,1],其中$x=0$和$x=1$为固壁边界。输出时刻为
  $t = 0.010$, $0.016, 0.026, 0.028, 0.030, 0.032, 0.034, 0.038$
\item[问题4]  
  \begin{displaymath}
    (\rho, u, p)(x,0)=\left\{
    \begin{array}{ll}
      (3.857143, 2.629369, 10.33333),& x <-4,\\
      (1 + 0.2 sin(5x), 0, 1),& x \geq-4,\\
    \end{array}
    \right.
  \end{displaymath}
  计算区间为[-5,5],其中在$x=\pm 5$边界处$\partial_x\rho=\partial_xu=\partial_xp=0$
  输出时刻为$t=1.8$
\end{description}
根据$E=\frac{1}{\gamma-1}p+\frac{1}{2}\rho u^2,\gamma=1.4$问题化为。
\begin{description}
\item[问题3]初始条件：
  \begin{displaymath}
    (\rho, \rho u, E)(x,0)=\left\{
    \begin{array}{ll}
      (1, 0, 2.5\times 10^3),&0\leq x <0.1,\\
      (1, 0, 2.5\times 10^{-2}),&0.1\leq x <0.9,\\
      (1, 0, 2.5\times 10^2),&0.9\leq x <1,\\
    \end{array}
    \right.
  \end{displaymath}
  计算区间为[0,1],其中$x=0$和$x=1$为固壁边界。输出时刻为
  $t = 0.010$, $0.016, 0.026, 0.028, 0.030, 0.032, 0.034, 0.038$
\item[问题4]  
  \begin{displaymath}
    (\rho, \rho u, E)(x,0)=\left\{
    \begin{array}{ll}
      (3.857143, 10.14185, 39.16666),& x <-4,\\
      (1 + 0.2 sin(5x), 0, 2.5),& x \geq4,\\
    \end{array}
    \right.
  \end{displaymath}
  计算区间为[-5,5],其中在$x=\pm 5$边界处$\partial_x\rho=\partial_x(\rho u)=\partial_xE=0$
  输出时刻为$t=1.8$
\end{description}



\section{数值方法}
\begin{description}
  \item[Lax-Friedrichs 格式($\lambda:=\frac{\tau}{h}$)]
    \begin{displaymath}
      u_j^{n+1}=\frac{u_{j+1}^{n}+u_{j-1}^{n}}{2}-\frac{\lambda}{2}(f(u_{j+1}^n)-f(u_{j-1}^n))
    \end{displaymath}
  \item[MacCormack 格式]
    \begin{displaymath}
      \left\{
      \begin{array}{rl}
        \bar{u}_j^*&=u_j^n-\frac{\tau}{h}(f(u_j^n)-f(u_{j-1}^n)).\\
        u_j^{n+1}&=\frac{1}{2}(u_j^n+\bar{u}_j^*)-\frac{\tau}{2h}(f(\bar{u}_{j+1}^*)-f(\bar{u}_{j}^*)).\\
      \end{array}
      \right.
    \end{displaymath}
\end{description}
\begin{description}
  \item[注1：]对于固壁边界条件中由于未给定边界上$E$和$\rho$的值,姑且假设他们保持常数，即与初值恒等。
  \item[注2：]对于问题4中的边界条件，迭代步中不计算边界点的值，一步迭代完成后，直接将边界点的值取为相邻点的值。
也就是用一阶差商代替微商。
\end{description}


\section{数值实验}

\subsection{LaxFriedrichs格式}
对问题中的$max\frac{\partial F}{\partial x}$的估计用了一个简单的试探的办法。
最后视$max\frac{\partial F}{\partial x}$为100，即取$\frac{\tau}{h}=0.01$以满足CFL条件。
LaxFriedrichs格式对间断的保持不好。
\begin{description}
  \item[问题3] 在程序中每计算到输出时刻就将$\rho$结果的值用plot函数作图输出。当网格不够小的时候，结果很不好。最后将网格加到20000个的时候得到下一章中的比较能接受的结果。
  \item[问题4] 同样当网格不够小时结果很难看。最后将网格加到了10000。
\end{description}

\subsection{MacCormach格式}
MacCormach格式没能计算出稳定的结果。在迭代几次之后，在间断附近出现强烈的震荡，然后震荡逐步扩散到整个区域。


\section{实验结果}
没有精确解，也没有合适的格式来解一个比较好的近似解，所以对计算结果只能有一个粗略的评价，没有办法做收敛阶的估计。
\subsection{LaxFriedrichs格式}
\begin{description}
  \item[问题3] 当网格数量在200的时候，计算出来的结果非常粗略。参照了\cite{thz03}的结果后，
    发现计算结果与精确解想去甚远，于是将网格大小加到20000得到如图\ref{im1}图\ref{im2}图\ref{im3}图\ref{im4}
    图\ref{im5}图\ref{im6}图\ref{im7}图\ref{im8}
    的结果
  \item[问题4] 这里给出网格数为2000（图\ref{im9}）和10000(图\ref{im10})的两个结果。可以看出网格数为2000的时候根本没有
    显示上面的波浪的形状。
\end{description}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_01.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.01时刻的计算结果}\label{im1}
\end{figure}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_016.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.016时刻的计算结果}\label{im2}
\end{figure}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_026.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.026时刻的计算结果}\label{im3}
\end{figure}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_028.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.028时刻的计算结果}\label{im4}
\end{figure}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_03.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.03时刻的计算结果}\label{im5}
\end{figure}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_032.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.032时刻的计算结果}\label{im6}
\end{figure}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_034.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.034时刻的计算结果}\label{im7}
\end{figure}

\begin{figure}[h!]\centering
\includegraphics[width=14cm]{./LF/q3_time_0_038.eps}
\caption{用LaxFriedrichs计算问题3在20000个网格的0.038时刻的计算结果}\label{im8}
\end{figure}

\begin{figure}[h!]\centering
  \includegraphics[width=14cm]{./LF/q4_2000.eps}
  \caption{用LaxFriedrichs计算问题4在2000个网格的1.8时刻的计算结果}\label{im9}
\end{figure}

\begin{figure}[h!]\centering
  \includegraphics[width=14cm]{./LF/q4.eps}
  \caption{用LaxFriedrichs计算问题4在10000个网格的1.8时刻的计算结果}\label{im10}
\end{figure}

\subsection{MacCormach格式}
MacCormach格式没能计算出稳定的结果。

\subsection{结论}
对于强间断的情况，单纯的守恒型格式经常得不到满意的结果。
\begin{thebibliography}{99}
% adaptive
\bibitem{thz01}
李治平，《偏微分方程数值解讲义》，2009
\bibitem{thz02}
汤华中，《计算流体力学讲义》，2011
\bibitem{thz03}
P. Woodward and P. Collella, J. Comput. Phys., 54(1984), 115-173.

\end{thebibliography}

\end{document}
